Energy of N Cooper pair by analytically solving Richardson- Gaudin equations 
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This Letter provides the solution to a yet unsolved basic problem of Solid State Physics: the 
ground state energy of an arbitrary number of Cooper pairs interacting via the Bardeen-Cooper- 
Schrieffer potential. We here break a 50 year old math problem by analytically solving Richardson- 
Gaudin equations which give the exact energy of these N pairs via N parameters coupled through 
N non-linear equations. Our result fully supports the standard BCS result obtained for a pair 
number equal to half the number of states feeling the potential. More importantly, it shows that 
the interaction part of the iV-pair energy depends on N as N(N — 1) only from N = 1 to the dense 
regime, a result which evidences that Cooper pairs interact via Pauli blocking only. 
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Superconductivity [D-Q is one of the most fascinat- 
ing phenomena of Solid State Physics. Its physical un- 
derstanding stayed a major problem for half a century 
[B| . The first step towards this understanding came from 
H. Frohlich Q who pointed out that two electrons with 
opposite spins can attract each other via the ion motion. 
A few years later, L. Cooper showed Q that, no matter 
how weak this attraction is, two electrons with opposite 
spins can form a bound state when added to a Fermi 
sea in an energy layer where a small attracting poten- 
tial acts - we will call this layer "potential layer" in the 
following. To handle more than one pair is clearly dif- 
ficult because this imposes to control Pauli blocking be- 
tween paired fermions. A way to overcome this difficulty 
is to turn to the grand canonical ensemble, with a pair 
number not fixed, as proposed by Bardeen, Cooper and 
Schrieffer (BCS). Using a wave function ansatz - based 
on the idea that electron pairs are bosonic particles, so 
that they are likely to condense all into the same state - 
they derived 0] , through a minimization of the Hamilto- 
nian mean-value, the condensation energy for electrons 
filling half the potential layer. 

A few years after BCS work, R.W. Richardson 0- 11 
and M. Gaudin 12, 13|, made an important step in show- 
ing that the BCS Hamiltonian leads to one of the very 
few exactly solvable Schrodinger equations. Indeed, the 
exact energy of N up/ down spin electrons paired by the 
so-called "reduced BCS potential", 
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with w p = 1 for e p in the potential layer (ep ,Ep + 0), 
reads as a sum of N complex quantities Rj solution of 
N algebraic equations. Although this is a significant 
advance compared to solving a second order differen- 
tial equation for ./V-body wave function, the resolution 
of these N equations still is a formidable math prob- 
lem which, over the past half century, stayed unsolved 
for arbitrary N and potential strength. Through an el- 



egant electrostatic analogy [3l , Richardson succeeded to 
recover the BCS condensation energy in the large N limit. 
For N small enough, these equations are commonly ap- 
proached numerically to understand the physics of super- 
conductor granules 15Ptl7 1 . 

Cooper pairs are composite bosons quite different from 
the semiconductor excitons we extensively studied over 
the last decade [Til] . Through our exciton studies, we 
however understood that the many-body physics of com- 
posite bosons is mainly driven by the Pauli exclusion 
principle between the particle fermionic components. In 
order to microscopically control the effect of Pauli block- 
ing when the Cooper pair number increases, it is nec- 
essary to stay in the canonical ensemble with both, the 
number of pairs and the number of states available for 
pairing, fixed. This is why we decided to tackle these 
Richardson- Gaudin equations again in order to solve 
them analytically for arbitrary TV. 

To grasp the trend induced by Pauli blocking on the 
pair binding energy, we first considered two pairs: even 
this N = 2 problem had no known solution although 
this definitely is the next problem to tackle after N = 1 
studied by Cooper. A year ago, we showed [l!| that the 
exact energy for two pairs reads as 
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Ei = 2e Fo —e c is the single pair energy found by Cooper. 
£p is the Fermi energy of the electrons which do not 
feel the potential. 2er/(l — cr) = £ c /f2 is the single pair 
binding energy in unit of the potential extension O. The 
dimensionless parameter a = exp(— 2/pV), with p being 
the density of states taken as constant in the potential 
layer, is sample volume free, p and 1/V linearly increasing 
with volume. 
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This 2-pair energy gives hints to understand the ef- 
fects of Pauli blocking on Cooper pairs. In a large sam- 
ple, p — > oo, so 9 — > 0. Difference between the ener- 
gies of two correlated pairs E 2 and two single pairs 2E± 
has a naive contribution 2/p which comes from the fact 
that the second pair must occupy the up/down states 
{ef + 1/p) just above the Ef Fermi level. E 2 — 2E\ also 
has a more subtle contribution 2e c /Na where Nq = pQ 
is the number of states in the potential layer from which 
paired electrons are formed. This 2e c /Nq contribution 
brings the 2-pair correlation energy from the two single- 
pair value 2e c down to 2e c {l — 1/Nq). This binding en- 
ergy decrease is induced by the moth-eaten effect com- 
ing from the Pauli exclusion principle between compos- 
ite bosons [lU . It is better understood by writting e c as 
Nqe* where e* = 2a/(l — a)p is the contribution of each 
of the Nq empty pair states in the potential layer. When 
a second pair is added, the number of states available 
for pairing decreases from Nq to Nq — 1; so the binding 
energy must decreases, due to Pauli blocking, from Nqe* 
to (Nq — l)e*, as we find. 

The purpose of this Letter is to see how this micro- 
scopic understanding of the Cooper pair correlation en- 
ergy extends to N > 2. We here present an analytical so- 
lution of Richardson-Gaudin equations for the energy of 
an arbitrary number N of up/down spin electrons paired 
by the BCS potential. We find that the iV-pair energy 
takes a remarkably compact form 



E N = NEi 



N(N - 1) 1 + a 



P 



1-(T 



(4) 



within terms in (N / p) n free from sample volume and thus 
negligible in the thermodynamic limit, in front of the two 
volume linear terms of En- 

This result reduces not only to Eq. for N = 2 but 
also to the energy obtained by Bardeen, Cooper and 
Schricffcr for Nbcs = Nq/2 which corresponds to a 
potential extending symmetrically on both sides of the 
normal electron Fermi sea. Indeed, Eq. (|4]) gives the con- 
densation energy as 
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where E 
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2 [Ne Fo + N(N - l)/2p] is the energy of 
A" "normal" electrons added above £f (1 , for a constant 
density of states p. This condensation energy reduces to 
EcNbcs/2 for half filling which is the BCS value, pA 2 /2, 
the excitation gap reading as A ~ Via 1 / 2 for a small. 

Eq.(4) is very astonishing at first because it looks as 
the first two terms of the small N expansion of the iV-pair 
energy: Higher order N terms seem missing! Actually, we 
have already reached these two first terms in a previous 
work j2(| in which we analytically solved Richardson- 
Gaudin equations in the very dilute limit N « N c where 
N c = pe c is the pair number over which single Cooper 



pairs would start to overlap. Since this number is far 
smaller than the BCS pair number Nq/2, the procedure 
we have used to extract the Rj 's from these equations is 
definitly not valid in the dense BCS limit. It moreover 
is so demanding that there were no hope to use it for 
higher order N terms in order to at least check that they 
do cancel exactly. 

To prove that these higher order terms do not exist 
in the large volume limit, we have constructed a totally 
different procedure. We have found a way to reach the 
sum of -Rj's directly, without calculating these parame- 
ters separately, as we have done in our previous work. 
This is somehow necessary because the Rj's are 2 by 2 
complex conjugate - with one Rj real for N odd. When 
A" increases, these Rj's run away from the real axis. This 
prevents their simple convergence for A~ large. By con- 
trast, J2j(Rj ~ Ei)/Nr] where 
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reduces to a degree-one rj polynomial due to a set of fun 



damental cancellations somewhat magic at first [21 



Since the analytical resolution of Richardson-Gaudin 
equations for arbitrary A~ and V is a formidable math 
problem which stayed open for 50 years, some readers 
may wish to see the major steps of our procedure. These 
are given below, in a separate section, to be easily skipped 
by a more general audience ready to accept that we have 
indeed proved Eq.(4). 

This very nice solution for a model hamiltonian widely 
used in the literature for BCS superconductivity but also 
cold gases, is not only missing in all textbooks but it 
mostly brings some interesting new light in a field com- 
monly considered as fully understood. Indeed, it provides 
a direct link between condensation energy and the Pauli 
exclusion principle, link out of reach when using the usual 
BCS ansatz in the grand canonical ensemble because N 
is not a free parameter in this ensemble. By rewritting 
the condensation energy for A^ pairs given in Eq.(5) as 

E^ — En = Ne c (N), we find that the average con- 
densation energy per pair in the Af-pair configuration, is 
proportional to the number of states available for pairing 
when (N — 1) pair states are already occupied. Indeed, 
for a total number of pairs Nq in this layer, the average 
condensation energy simply reads 



£ C (N) = [Nq — (AT — 1)] e*. 
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This result evidences that the unique consequence of in- 
creasing the number of pairs in the potential layer is to 
block more and more states in this layer, until all states 
are occupied: the correlation energy would then reduce to 
zero. In the Cooper problem, N = 1 and e c (l) = Nqe* 
while in the usual BCS configuration A^ = Nq/2 and 
the average condensation energy reduces to one half the 
single pair value. Eq.(7) shows that Cooper pairs only 
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"interact" through Pauli blocking. As a direct conse- 
quence, they can overlap without breaking, by contrast 
to excitons which dissociate into an electron-hole plasma 
through a Mott transition when overlap starts. 

One important consequence of Eq.(7) is that the aver- 
age pair energy cannot be identified, as commonly done 
in textbooks [2| , with the excitation gap much larger than 
the single pair value. The gap is the energy to break a 
pair. In the dense regime, this energy is far larger than 
just the energy of the broken pair due to Pauli induced 
many-body effects between the broken pair and the re- 
maining (N — 1) unbroken pairs. 

Mathematical resolution of Richardson-Gaudin 
equations. The exact energy of N up and down 
spin electrons interacting via the BCS potential given 
in Eq.(l), reads 0, [TJ, \2% as £)" 
are coupled through 
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We now list the major steps of our resolution. 

• The first step is to note that for one pair, E\ = 
R\, fulfills the above equation without the k sum. By 
substracting this equation from Eq. (8), we get 
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The potential V has disappeared from this equation: it 
is now hidden into E\ . We then set Rj — E± = e c Zj and 
expand the sum over p as a Zj infinite series. By setting 
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we find, for 7 defined in Eq. (|3j), that the Zj's fulfill 
00 .. 
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The unpleasant part of this equation is its second term. 

• To get rid of it, we, in a second step, multiply the 
above equation by z l j with I = (0,1,2,...) and we sum 
over j. This gives 



= ^2 a n Z n+ i + jDi, 



where the sums Zm and D m are defined as 
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The iV-pair energy follows from Z\ . 

Through the replacement of zj 1 by (z™ — z™)/2 in 
D m: it is easy to show that D Q = 0, D x = N(N -l)/2, 



D 2 = (N-1)Z U D 3 = (N-3/2)Z 2 + Zf/2, and so on... 
As Zq = N, the general expression of D m actually reads 
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• The third step is to rescale Z m as Z m = Nrj m X m 
with r\ defined in Eq.(6). Eqs. (12-14) then give the set of 
equations fulfilled by the X m 's as 
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for £ ^ 1 while for £ = the RHS reduces to zero. This 
readily shows that the X m 's are rj series. Let us write 
them as X m = J2 q ^o x m, q Tf- Since X = 1, we get 
xo.o = 1 and xo, q =£o = 0. The other x m> qS follow from 
identification of the rf terms in Eq.(15). A tedious but 
straightforward calculation shows that X\ is an odd func- 
tion of r\ 

where X[ depends on (77, 1/N,a) as 
X[ = {l-^)(l+a)(l-*Y{yW + (y 5 , V -f)v 5 

■}, (17) 

the j/m.fc's depending on a only. This leads to 
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Since the last term scales as (rf ,77 s ,rf , • • ■), it gives 
volume free contributions to En in (N / p) n with n — 
(2, 3, ■ ■ • ), in agreement with Eq.(4). 

• To be complete, the last important step of our res- 
olution is to prove that all corrections to the first two 
terms of E^ are indeed volume free, i.e., X\ reduces to 



(0) 



-77(1 + c)/4 for large volume. To do it, we re- 



consider Eq.(15) without its last term 



ai X^ + a 2r] xi% + a 3 ifx(% + ■ ■ ■ = \ £ Ifx|l r 



(19) 

The ll°''s, solution of this equation, are rj polynomials 

-y-(O) _ T (o) , T (o) , , (0) m (2()) 

^-m — x m,0 ' x m,l'l ' ' x m,m'l • I 

the x^ k coefficients, defined for ^ k ^ m and nonzero 
for even m + k only, being just the ones appearing in 

P ( ta ) = " u t _ (T \ni t -l) n ] (21) 

y ' 2 n (n+ 1)! dt n LV ' y ' 1 v ' 

- T (o) , T (o) f+ , (0) , n (22 ) 
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From Pi(t,a) = -(1 + <r)/4 + t/2, we then find x[ 0) = 

x^\rj = —77(1 + cr)/4 as we want. 

The first key to check that the solution of Eq.(19) is 
given by Eq.(20) is to note that P n (t,a)'s is related to 
the Legendre polynomials L n (t) through 



P n {t,a) = 
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The orthogonality of Legendre polynomials |23j then 
gives for $C k ^ n — 1 



1 r 1 

= / t k P n (t,a)dt. 



(24) 



By writting P n (t,a) as in Eq.(22), it follows from this 
equation that = x$ n a k+ i + x| l °| lj „_ 1 a fe+2 + • • • + 



c (0) 

L 2n,0 



ifc+i+n, where a n is just the scalar appearing in 
Eq.(lO). 

The second key also follows from the link between 
P n (t, it)'s and the L n (t) Legendre polynomials. It is pos- 
sible to show that 
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the RHS reducing to zero for n = 0. Eq.(22) inserted 
into the LHS of Eq.(25) gives this LHS as a t poly- 



nomial G n (t) = a;„+i i „_ifl , i(*) + 1" x 2n,o9n(t) where 

g m (t) again depends on the aj scalar defined in Eq.(10) 
as g m (t) = Y^,j=i ajt m ~ J ■ For n = I, this in particular 

gives G\(t) = aix^Q = 1/2 since Po{t,a) = 1. 

Using these two equalities, it becomes possible to show 
by identification that Eq.(20) fulfills Eq.(19). Details will 
be given in an extended version of this Letter. 

Conclusion: We here derive the energy of N Cooper 
pairs by solving Richardson-Gaudin equations analyti- 
cally for arbitrary A and potential strength. We prove 
that, for large samples, the interaction part of the A- 
Cooper pair energy depends on A as A(A — 1) only: 
higher order A terms do not exist. As a result, the aver- 
age Cooper pair binding energy (E$ — E^)/N linearly 
decreases with pair number from N = 1 to the dense 
regime, (sec Eq.(5)), this energy being simply propor- 
tional to the number of empty states available for pair- 
ing in the potential layer. As a result, the pair binding 
energy cannot be identified with the gap, as commonly 
done. 

Our result fully supports the BCS result for the ground 
state energy obtained in the grand canonical ensemble 
when the potential layer extends symmetrically on both 
sides of the normal electron Fermi sea. One puzzling 
question still remains: Why the BCS ansatz leads to the 



exact ground state energy for N equal to half-filling since 
its projection onto the iV-pair subspace corresponds to 
\ipj^ C }=(B^) N \F ) with all pairs condensed into the 
same state, while the exact Richardson-Gaudin ground 
state reads as \tp N ) = B^(Ri) ■ ■ ■ (R N )\F ) where 
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the i?j's being all different due to Pauli blocking, as 
seen from the last term of Eq.(8)? Is this A-pair energy 
agreement also valid for correlation functions? We hope 
that this Letter will stimulate more works in connexion 
with the effects of the Pauli exclusion principle on paired 
electrons, in a field, conventional BCS superconductivity, 
commonly considered as fully understood. 

We wish to thank Tony Leggett, Walter Pogosov, and 
Guojun Zhu for numerous discussions. 
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